% this funciton will read the data stored in meta image files (*.MHA)
function im = mharead(varargin)

if nargin < 1
    % user did not specify an path, give a dialog to do so
    % this will store the last path name in the main app data structure
    lastPath = getappdata(0, 'lastPath');
    if isempty(lastPath); lastPath = 'C:\'; end;
    [name, pathstr] = uigetfile('*.mha', 'Pick meta image to open', lastPath);
    setappdata(0, 'lastPath', pathstr);
    mhapath = fullfile(pathstr,name);
elseif nargin == 1
    % user specified a path to an mha file
    % break down input file name
    mhapath = varargin{1};
    [pathstr, name, ext, versn] = fileparts(mhapath);
elseif nargin > 1
    % i do not know what to do with this many parameters
    disp('ERROR: This function expects 0 or 1 arguments.')
    return
end

fid  = fopen(mhapath);
mhatext = textscan(fid,'%s','delimiter','\n');
fclose(fid);

% parse headers, assumes that tags are always in the same order.
% do i need to check the tag to make sure??

for t = 1:length(mhatext{1})
    tag = textscan(mhatext{1}{t},'%s',1);
    tag = char(tag{1}{1});
    switch lower(tag)
        case 'ndims'
            nDims = sscanf(mhatext{1}{t},'%*s%*s%f');
        case 'elementspacing'
            elementSpacing = sscanf(mhatext{1}{t},'%*s%*s%f%f%f');
        case 'dimsize'
            dimSize = sscanf(mhatext{1}{t},'%*s%*s%f%f%f');
        case 'elementtype'
            elementType = sscanf(mhatext{1}{t},'%*s%*s%s');
        case 'elementbyteordermsb'
            elementByteOrderMSB = sscanf(mhatext{1}{t},'%*s%*s%s');
        case 'elementdatafile'
            elementDataFile = sscanf(mhatext{1}{t},'%*s%*s%s');
        case 'anatomicalorientation'
            AnatomicalOrientation = sscanf(mhatext{1}{t},'%*s%*s%s');
        case 'centerofrotation'
            CenterOfRotation = sscanf(mhatext{1}{t},'%*s%*s%f%f%f');
        case 'offset'
            Offset = sscanf(mhatext{1}{t},'%*s%*s%f%f%f');
        case 'transformmatrix'
            TransformMatrix = sscanf(mhatext{1}{t},'%*s%*s%f%f%f%f%f%f%f%f%f');
        case 'binarydata'
            BinaryData = sscanf(mhatext{1}{t},'%*s%*s%s');
        case 'objecttype'
            ObjectType = sscanf(mhatext{1}{t},'%*s%*s%s');
        otherwise
            disp(['MHA field not recognized: ' tag])
    end        
            
end


% need to handle 2 and 3d images
if length(dimSize) == 2; dimSize(3) = 1; end
% initialize image size
im(dimSize(2),dimSize(1),dimSize(3)) = 0;

% need to convert element type to matlab precision format
switch elementType
    case 'MET_CHAR'
        elementType = 'schar';
    case 'MET_UCHAR'
        elementType = 'uchar';
    case 'MET_SHORT'
        elementType = 'short';
    case 'MET_USHORT'
        elementType = 'ushort';
    otherwise
        disp('ERROR: Unknown element type.')
        return
end
        

% now that we have the particulars, we can read the data file
if ~strcmp(elementDataFile,'local')         % Determine whether the binary data is local or stored in another file
    bipath = fullfile(pathstr,elementDataFile);
    fid = fopen(bipath);
    [fname, mode, mformat] = fopen(fid);
    for j = 1:dimSize(3)
        for i = 1:dimSize(2)
            im(i,:,j) = fread(fid,dimSize(1),elementType,mformat);
        end
    end
end
fclose(fid);    
return
